#########################################################
# SIMULATIONS OF ALTERNATIVE DENOMINATOR SPECIFICATIONS #
#########################################################

# Author: Kasia Nalewajko
# First version: 29 November 2023
# Replicated: 13 June 2024


rm(list = ls())

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/jpop_simulations.Rda")
load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/main.Rda")

# RUN MODELS --------------------------------------------------------------

# Run 1,500 regressions recording p-values of the main independent variable
# Plot the p-values

simulations <- left_join(main, for_simulations)
temp <- simulations %>% 
  mutate_at(vars(starts_with("jpop_random_distribution_")), funs(logged_random_proportion = log((all_jewish_victims+1)/(.+1)+1)))

vars_with_random_distributions <- for_simulations %>% 
  dplyr::select(starts_with("jpop_random_distribution_"))

# Create empty containers
dvars <- grep(x = names(temp), pattern = "_logged_random_proportion$", value = TRUE)
p <- rep(NA, 1500)
models_list <- list()

# Loop over outcomes
for (i in 1:length(dvars)) { # denominators
  
  models_list[[i]] <-  fixest::feols(as.formula(paste(dvars[i], "~ + log(", vars_with_random_distributions[i], "+1) + log(pop1936+1) + synagogues_dummy + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1)  + catholic_churches + log(FRANCISTE_36_1_pct+1) + log(ActionFrancaise_19_pct+1) + turnout_36_1_pct + log(nuance_droite_36_1_pct+1) + log(nuance_centredroit_36_1_pct+1) + log(nuance_centregauche_36_1_pct+1) + log(nuance_gauche_36_1_pct+1) + log(nuance_extremegauche_36_1_pct+1) + area_sqkm + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name) | log((FFIFFCgentile*1000/pop1936)+1) ~ log(mpf1000+1)")),
                                     cluster = ~c(id_bureau, resregion2), 
                                     data = temp)
  
  dvars[[i]] <- dvars[i]
  p[[i]] <- models_list[[i]][["coeftable"]][1,4]
  
  print(i)
}

pvalues_simulatedprops <- data.frame(
  dvars = dvars,
  p = p)

# PLOT -------------------------

pvalues_simulatedprops %>% 
  ggplot() +
  aes(x = p) %>% 
  geom_histogram(binwidth = 0.01) +
  geom_vline(aes(xintercept = 0.05), linetype = 1, color = "gray20") +
  geom_text(aes(x=0.05, label="0.05\n", y=100), colour="blue", angle=90, text=element_text(size=11)) +
  theme_bw() +
  labs(
    x = "P-values",
    y = "Count regressions"
  )

# EXPORT ---------

ggsave(
  "D1_pvalues_simulatedprops.png",
  plot = last_plot(),
  path = "./00 SUBMITTED/00 APSR final/03 dataverse_online_appendix/figures",
  width = 21,
  height = 18,
  units = "cm",
  dpi = 300
)

